# misc.funE_B -- 
#
# Copyright (c) 1988,2001, Willem Jan Zaadnoordijk

from math import exp,log,cos,sqrt

def E1(u):
	"exponential integral"
	a1=-.57721566
	a2=.99999193
	a3=-.24991055
	a4=.05519968
	a5=-.00976004
	a6=.00107857
	af1=.250621
	af2=2.334733
	bf1=1.681534
	bf2=3.330657
	if (u <= 1) :
		rfe1=((((a6*u+a5)*u+a4)*u+a3)*u+a2)*u+a1
		rfe1=rfe1-log(u)
	elif (u <= 111.):
#     prevent overflow in DEXP (111. large enough not to lose accuracy
#     after the conversion to single precesion: DEXP(111)=1.609487e+48)
		rfe1=((u+af2)*u+af1)/((u+bf2)*u+bf1)
		rfe1=rfe1/(exp(u)*u)
	else :
		rfe1=0.
	return rfe1


def erfc(x) :
	"complementary error function erfc(x)= 1-erf(x)"

#  using formula 7.1.26 on page 299 of Abramowitz & Stegun
	p = 0.3275911
	a1= 0.254829592
	a2=-0.284496736
	a3= 1.421413741
	a4=-1.453152027
	a5= 1.061405429

	xabs=abs(x)

	t=1./(1.+p*xabs)
	rferfc=t*(a1+t*(a2+t*(a3+t*(a4+t*a5))))*exp(-x*x)
	if (x < 0.):
		rferfc=2.-rferfc

	return rferfc


def J0(x) :
	"Bessel function first kind order zero"
	a0= 1.
	a1=-2.2499997
	a2= 1.2656208
	a3=-.3163866
	a4= .0444479
	a5=-.0039444
	a6= .00021
	b0= .79788456
	b1=-.00000077
	b2=-.0055274
	b3=-.00009512
	b4= .00137237
	b5=-.00072805
	b6= .00014476
	c0=-.78539816
	c1=-.04166397
	c2=-.00003954
	c3=.00262573
	c4=-.00054125
	c5=-.00029333
	c6=.00013558

	if (abs(x) <= 3) :
#     Abramowitz & Stegun (1972) page 369, (9.4.1)
		r=x*x/9.
		rfj0=a0+r*(a1+r*(a2+r*(a3+r*(a4+r*(a5+r*a6)))))
	else :
#     Abramowitz & Stegun (1972) page 369-370, (9.4.3)
		r=3./x
		f=b0+r*(b1+r*(b2+r*(b3+r*(b4+r*(b5+r*b6)))))
		t=x+c0+r*(c1+r*(c2+r*(c3+r*(c4+r*(c5+r*c6)))))
		rfj0=f*cos(t)/sqrt(x)
	return rfj0


def J1(x) :
#     Bessel function first kind order one
	a0=  .5
	a1= -.56249985
	a2=  .21093573
	a3= -.03954289
	a4=  .00443319
	a5= -.00031761
	a6=  .00001109
	b0=  .79788456
	b1=  .00000156
	b2=  .01659667
	b3=  .00017105
	b4= -.00249511
	b5=  .00113653
	b6= -.00020033
	c0=-2.35619449
	c1=  .12499612
	c2=  .0000565
	c3= -.00637879
	c4=  .00074348
	c5=  .00079824
	c6= -.00029166
	if (abs(x) <= 3) :
#     Abramowitz & Stegun (1972) page 370, (9.4.4)
		r=x*x/9.
		rfj1=(a0+r*(a1+r*(a2+r*(a3+r*(a4+r*(a5+r*a6))))))*x
	else :
#     Abramowitz & Stegun (1972) page 370, (9.4.6)
		r=3./x
		f=b0+r*(b1+r*(b2+r*(b3+r*(b4+r*(b5+r*b6)))))
		t=x+c0+r*(c1+r*(c2+r*(c3+r*(c4+r*(c5+r*c6)))))
		rfj1=f*cos(t)/sqrt(x)
	return rfj1


def J00i(i) :
	"function returns Ith zero of J0"
#  first twenty from Abramowitz & Stegun (1972) table 9.5 page 411
#  second twenty from Watson (1958) table VII page 748
	zero=[
		 2.4048255577,  5.5200781103,  8.6537279129, 11.7915344391,  \
		14.9309177086, 18.0710639679, 21.2116366299, 24.3524715308,  \
		27.493479132 , 30.6346064684, 33.7758202136, 36.9170983537,  \
		40.0584257646, 43.1997917132, 46.3411883717, 49.4826098974,  \
		52.6240518411, 55.765510755 , 58.9069839261, 62.0484691902,  \
		65.1899648, 68.3314693, 71.4729816,  74.6145006, 77.7560256, \
		80.8975559, 84.0390908, 87.1806298,  90.3221726, 93.4637188, \
		96.6052680, 99.7468199, 102.8883743,106.0299309,109.1714896, \
		112.3130503,115.4546127,118.5961766,121.7377421,124.8793089 ]

	if (i > 40) :
		print 'error: function J00i, only 40 zeros available, i too large'
		return
	else :
		return zero[i-1]
